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Abstract 

Various Monte Carlo programs, developed either by small groups or widely available, have been 
used to calculate the effects of decays of radioactive chains, from the original parent nucleus to 
the final stable isotopes. These chains include uranium, thorium, radon, and others, and generally 
have long-lived parent nuclei. Generating decays within these chains requires a certain amount 
of computing overhead related to simulating unnecessary decays, time-ordering the final results in 
post-processing, or both. We present a combination analytic/stochastic algorithm for creating a 
time-ordered set of decays with position and time correlations, and starting with an arbitrary source 
age. Thus the simulation costs are greatly reduced, while at the same time avoiding chronological 
post-processing. We discuss optimization methods within the approach to minimize calculation 
time. 
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i I. INTRODUCTION 



2 Various classes of physics experiments require reliable simulations of radioactive sources. 

3 Obvious examples include individual radiation detectors, but less obvious examples are low- 

4 background experiments such as are used to search for neutrinoless double-beta decay or 

5 WIMP Dark Matter. The radioactive sources are predominantly either natural background 
e radiation that may obscure any potential signal, or calibration sources used to experimen- 

7 tally determine detector performance. Given the ultra-sensitive nature of current and next- 

8 generation experiments, previous understanding of detector response to bulk material ra- 

9 dioactivity and calibration sources is usually not sufficient to meet stated scientific goals. 

10 Simulations of these detectors therefore require event generators that are highly accurate in 

11 temporal and spatial extent. 

12 Various Monte Carlo programs are available to calculate the effects of radioactive decays 

13 in a detector setup. These programs include GEANT4 [1], MCNP [2], FLUKA [3], or 

14 perhaps other programs written by smaller groups or individuals. Often, the radioactive 

15 decays in question derive from the decay of a single isotope to a stable nucleus, such as 

16 in the cases of 137 Cs, 60 Co, and 40 K. In other cases, however, a chain of decays may occur 
l? before a stable nucleus is reached, as in the cases of 238 U, 232 Th, various isotopes of radon, 
is and so on. Within these decay chains, temporal and spatial correlations are often relevant. 

19 One reason is to determine the effects of detector dead time on the acquired signal. Another 

20 reason is to interpret multiple events within any given time frame and energy range as being 

21 correlated via a single decay chain, and therefore not contributing to the signal rate. 

22 In running Monte Carlo simulations, whether for natural backgrounds or calibration 

23 purposes, the age of the source needs to be taken into account. If the source is young and 

24 the chain has not reached secular equilibrium, the late decays will have less activity than if 

25 the chain is in secular equilibrium. In the case of background environmental radiation, the 

26 decay chains are often taken to be in secular equilibrium, and the effective source age can be 

27 on the order of millions of years. In the case of manufactured sources, the age will usually 

28 be much shorter, e.g. up to 20 years or so in the case of 228 Th. 

29 One simple and straightforward method to account for source age is to simulate many 

30 decays, starting with the top-most isotope, and record the time of each decay in the chain. 

31 Meanwhile, the time of each subsequent top-most decay is determined using that isotope's 
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1 activity. In the analysis step, a cut can be made rejecting all events before the source age. 

2 If this approach is used, the number of decays required to achieve secular equilibrium can 

3 be very large. For example, the longest-lived descendant of 238 U is 234 U, with a half-life of 

4 245.5 thousand years. Secular equilibrium is obtain after approximately 10 half-lives of the 

5 longest-lived daughter, which in this example means 2.5 million years of simulated events. If 
e the decay activity were taken to be a mere 1 mBq, this approach still requires on the order 
? of 10 13 decays to process — an untenably large number of decays, even by modern computer 
a standards. 

9 It may be possible to solve this issue by incorporating an analysis of detector response 

10 time, although there are a couple disadvantages to this approach. One is that any event 

11 generator created in this manner would be applicable only for that specific kind of detector. 

12 Another is that the total source activity must also be taken into account; when simulating 

13 low-background events, the event generators may incorporate one detector response time, 

14 while simulating calibration events may require a different detector response time. The event 
is generator would become even more specialized, and lose its general applicability. 

16 A second difficulty comes in post-processing the data to put it in chronological order. 

17 Continuing the above example of a 238 U source, the decay of one isotope in a source is not 
is necessarily directly followed in time by the decay of its own daughter. Rather, the next decay 

19 might derive from a different 238 U nucleus. For this reason, there may be effects of detector 

20 pileup or dead time that cannot be taken into account in the simulation by simply following 

21 the time structure of decay from a single parent to its eventual stable descendant. To solve 

22 this problem, the data can be post-processed so that all decays are placed in chronological 

23 order, regardless of their individual ancestor isotopes. A single simulation run can involve 

24 billions of energy depositions, and the computer resources necessary to sort these events 

25 may be substantial. 

26 Accurate temporal correlations are, of course, not always the top priority in a detector 

27 simulation. Historically, the background levels of various materials have not always been 

28 known, and the simulation data rate simply ends up being scaled to the experimental data 

29 to infer the physical source rates. In such situations, the time of the events is usually not 

30 incorporated into the analysis. Once the absolute levels of radioactivity are known, however, 

31 more accurate simulations can be run with the correct inter-event timing. 

32 In this paper, we present an algorithm to obtain time-ordered events from a decay chain 
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1 with an arbitrary source age, while preserving temporal and spatial correlations. Our ap- 

2 proach incorporates data structures to minimize the time required for both random search 

3 of the accumulated history of many decays, as well as insertion and deletion operations. 

4 We also present methods to reduce the recorded decay history data, minimizing computer 

5 memory usage. We evaluate our results using the 232 Th and 228 Th decay chains. We also 
e benchmark the time and memory usage for various starting activities and total number of 
? desired decays in the data set. 

8 II. COMBINING STOCHASTIC AND DETERMINISTIC METHODS FOR DE- 

9 CAY CHAIN GENERATION 

10 A purely stochastic approach to decay chain generation is intuitively compatible with 

11 Monte Carlo simulations, and its results have a high degree of reliability. In such an ap- 

12 proach, a top- most nucleus, e.g., 238 U for the 238 U chain, decays, followed by a daughter. 

13 The daughter decays after a time randomly sampled from the appropriate exponential time 

14 distribution with the well-known form 

Probability ~ e -* fa ( 2 )/ T i/2 (1) 

is where Ti/ 2 is the half life of the parent. The daughter itself then decays, with the time until 

16 the next decay determined by the next daughter's half life, and so on until a stable nucleus 

17 is reached. The time between top-most decays is sampled from the distribution e~ tr , where 
is r is the activity of decays of the top-most radionuclide. This approach, while conceptually 

19 simple, often cannot be implemented because of the concerns regarding age of the source 

20 and computing requirements discussed in the previous section. 

21 A second, more subtle problem with the purely stochastic approach is the limitations of 

22 standard data types. The double floating point type, for example, only has 16 decimal digits 

23 of accuracy. Over the course of the 2.5 million years it takes to obtain equilibrium in a 238 U 

24 decay chain, a double floating type can maintain accuracy only to approximately 10 ms. 

25 Yet the shortest-lived daughter in the chain, 214 Po, has a half-life of 164.3 /is. In a purely 

26 stochastic implementation, a 214 Po decay would occur simultaneously with its parent, 214 Bi, 

27 even though temporal separation of hundreds of microseconds is well within the capabilities 

28 of many physical detectors today. While this discrepancy between simulation and experiment 
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1 may seem trivial, resolving the issue would alleviate any concerns over this artificial pileup 

2 of events. 

3 It may be enticing to use a deterministic approach to calculate the decays from a chain. 

4 The equations involved have long been solved by H. Bateman [I]. Unfortunately, a purely 

5 analytic approach is incompatible with a Monte Carlo simulation of radioactivity, as the 
e individual decays lose their spatial correlations. Consider the case of a 212 Bi nucleus decaying 

7 into a 212 Po nucleus. Because 212 Po has a half life of 300 ns, its own decay will be tightly 

8 correlated to the position of the 212 Bi, and thus information from each decay must be passed 

9 along to progeny nuclei. This information is lost in the purely deterministic approach. 

10 We have therefore combined the two approaches in the current treatment. We employ 

11 Bateman's equations to calculate populations of all isotopes in the decay chain at any given 

12 source age. With those populations known, we proceed using a stochastic approach. In this 

13 way, we avoid the severe computational overhead while preserving the spatial correlations 

14 within a decay chain. 

15 A. The algorithm 

16 The inputs to the algorithm for any given decay chain, apart from defining the chain itself, 

17 are the starting activity of the top-most isotope, the age of the source, and the number of 
is events desired in the resulting decay record. Using these inputs makes direct comparison to 

19 known sources more straightforward, as commercial radioactive sources have an activity and 

20 manufacturing time stamp. In the case of environmental radiation, such as background 238 U 

21 or 232 Th, the steady-state activity is usually measured in situ, and the age of the source then 

22 set to 10 times the half life of the longest-lived daughter nucleus to reach secular equilibrium. 

23 The first step in the process is to calculate the starting populations of all isotopes in the 

24 chain using the Bateman equations. In our own implementation, we used the Moral and 

25 Pacheco approach to the Bateman equations [5j. Following Moral and Pacheco's recommen- 

26 dation, we used Mathematica [6] to solve the equations, and entered the solutions into our 

27 decay generator computer code. We assumed, as in [3], that the population of the chain 

28 before source aging was entirely in the top- most isotope. 

29 Because the population equations are analytic, it is possible to obtain a fractional pop- 

30 ulation. To limit calculations performed on unphysical numbers, if the population ever fell 

5 



1 below 0.5 in our treatment, we adjusted the population to identically 0. Once the popula- 

2 tions were determined, the starting activity of the i th isotope, r i; was calculated using the 

3 equation 



4 where iVj is the population of the i isotope at any given time. 

5 Within the current algorithm, we distinguished between decays from the primary pop- 
e ulation, i.e., the population after source aging, and decays from the subsequent chains. 
7 Additionally, we distinguished and tracked the time of primary decays ("primary time") 
e and the time of chain decays ("chain time"). The algorithm is as follows: 

9 1. For the calculated surviving primary populations, determine rates for each isotope 

10 using Eq. (J2J) 

11 2. Increment the primary time based on an exponential random sampling, using the total 

12 surviving primary decay activity in Eq. ([!]) 

13 3. Randomly select the next primary decay isotope, weighted by the activities of the 

14 surviving primary populations 

15 4. Reduce the primary population of that isotope by 1 

16 5. Set the chain time to the primary time 

17 6. Generate a location for the new primary decay 

is 7. Add the decay to the record history, including isotope, position, and time 

19 8. If the isotope has a radioactive daughter, increment the chain time by a random 

20 number distributed according to Eq. ([TJ, with the appropriate Ty 2 of the isotope, and 

21 propagate the decay location to the new decay 

22 9. If the number of decays in the data record is greater than the number requested, 

23 remove the last entry from the record 

24 10. Continue the previous three steps until a stable isotope is reached 

25 11. Return to the first step unless the requested number of primary decays has been 

26 reached 

27 As can be seen from this algorithm, the correct spatial and temporal correlations through 

28 an entire decay chain, starting at an arbitrary source age, is maintained. Additionally, the 
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1 issue of temporal resolution is largely resolved. Consider the case of accumulating 10 million 

2 decay events in the data record. If the total activity were a constant 1 Becquerel, then the 

3 last decays in the data set would have a time resolution on the order of 1 ns. 

4 Some care must be taken in determining the correct termination conditions. If 10 million 

5 events are desired in the decay record, it is incorrect to simply stop the calculation after 10 
e million events have been entered. 232 Th, for example, has 10 products in its decay chain, so 
? 1 million decays of 232 Th would provide the desired number of events. Unfortunately, after 
s those 1 million decays of 232 Th, with a relatively steady activity, the event rate will drop 

9 precipitously, as the only decays after that point in time depend on the 5.75-year half life of 

10 228 Ra. To avoid this unphysical data record, 10 million full decays of the starting population 

11 should be performed. 

12 The termination conditions can also be optimized. For example, if the number of re- 

13 quested decays has been entered into the data record, and the decay time is beyond the 

14 time of the last decay in the record, there is no point in continuing that particular chain, 
is Other optimizations may also exist, depending on the specific circumstances of starting 

16 activities, number of requested events, and so on. 

17 As a comment on the algorithm, it may be more natural when running a simulation to 
is require a set amount of simulated time, rather than a set number of simulated events. Unfor- 

19 tunately, this introduces a kind of tension with, for example, GEANT4, where the number 

20 of events is specified as the input. To run for a specific length of time, we recommend calcu- 

21 lating beforehand how many events there should be in the desired data set, and augmenting 

22 that number appropriately to ensure that the full desired simulation time is accumulated in 

23 the event record. 

24 B. Data structure 

25 Our approach required an appropriate data structure to maintain and time-order the 

26 events as they were generated by the algorithm. The standard C++ sequence container 

27 classes of vectors or lists are disadvantageous because of computing overhead. Vectors have 

28 the advantage of random access, allowing for a find operation with 0(log n) overhead. 

29 Unfortunately, insertion and deletion operations are computationally expensive, with O(n) 

30 requirements. Conversely, lists handle insertion and deletion in 0(1) time, but without 
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1 random access, finding the correct location to place a decay is an O(n) operation. We 

2 therefore decided to use a binary search tree data structure, where a find operation is 

3 performed in 0(log n) time, while insertion and deletion are still maintained at 0(1), thus 

4 combining the advantages of both vectors and lists. The C++ associative container class 

5 multimap would also be appropriate, although we implemented our own binary search tree 
e to avoid residual overhead of a multimap. 

7 Each subsequent primary decay occurs later in time than the one before it. Similarly, 

s a chain decay also occurs after its parent decay. These are obvious physical facts, but can 

9 negatively affect the performance of the data structure. Because each decay time occurs 

10 after the one before it, it is possible to end up with a semi-degenerate binary search tree, 

11 far more heavily weighted to one side than the other. If this were to occur, find operations 

12 approach O(n), and we lose the advantage of the search tree. 

13 To fix this situation, we pre-seeded the search tree with empty nodes to install an evenly- 
u spaced structure, and thus avoid a quasi-degenerate tree. We calculated the anticipated 

15 simulated time window to acquire all requested decays. This time window was approximate, 

16 as the ever-reducing primary population altered the time of the last decay in a non-trivial 

17 way. We then sub-divided this time window into a number of levels, where each level 
is contained twice as many empty nodes as the level above it. The spacing of each level was 

19 such that the node split in half the time to either side of the upper nodes (see Fig. [I]). The 

20 levels also had to be inserted in the order of their level, again to avoid quasi- degeneracy. 

21 The time savings of establishing an empty node structure in the binary search tree varied 

22 with the number of events already entered. If only one level was used, the time it took to 

23 populate the tree with 10 5 decays took approximately ten times as long as if 20 levels were 

24 used. As the tree is populated, however, it becomes increasingly balanced, since primary 

25 decays tend to occur at earlier times and chain decays at later times. This leads to a counter- 

26 intuitive result that the time savings afforded by pre-seeding the tree diminish with an 

27 increasing number of event. The computing time required to populate a tree with 10 7 events 

28 did not change appreciably with the number of pre-seeded levels. The authors therefore 

29 recommend using 20 pre-seed levels, or approximately 1 million empty nodes, to provide 

30 adequate protection against degeneracy, at the expense of less than 2 seconds of computing 

31 time and 40 megabytes of memory. If a large number of independent simulations, each with 

32 10 5 decays, is required, the time savings might be substantial. 
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FIG. 1. The structure of empty nodes in the pre- 
seeded binary search tree. By subdividing the time 
window, degeneracy in the data structure can be 
minimized. The nodes were entered in the order of 
their level number, and split the time between the 
two nodes closest in time, e.g., the time of the node 
circled in blue is at the average of its parent node 
and the left node of Level 2. 



i III. APPLICATION AND EVALUATION 



2 We evaluated our approach using the 232 Th and 228 Th decay chains. The former intro- 

3 duces a case of a near-constant decay activity in the top-most isotope, while in the latter 

4 case the population of the parent nucleus noticeably decays with time. 



5 A. Using a 232 Th source 

e Figure [2] show the decay activities of all the isotopes in the 232 Th decay chain. For 

7 the three figures, the source age was set to 0, 1, and 5 times the half life of the longest- 

a lived daughter in the chain, 228 Ra. The figures demonstrate the reliability of combining the 

9 analytic Bateman calculations to determine the starting population, followed by stochastic 

10 decays within the chain. For these calculations, the starting decay activity of the 232 Th was 

11 kept constant at 0.3 mBq, and 10 million decays were recorded. 

12 Figure [2] also shows how the current approach also properly handles branching within the 

13 decay chain. In this case, 212 Bi decays to 208 T1 (branching ratio = 35.93%), 212 Po (branching 
u ratio = 64.05%), and directly to 208 Pb (branching ratio = 0.023%). The rates of 208 T1 and 
15 212 Po are therefore lower than the others. 
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FIG. 2. Decay activities of isotopes within the 
232 Th decay chain. The top figure has a source 
age of years, the middle a source age of 1 half 
life of 228 Ra, and the bottom figure a source age 
of 5 half lives of 228 Ra. In each case, the vertical 
dashed line marks 10 half lives from the creation of 
the source. For reference, the half life of 228 Ra is 
5.75 years. 10 
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FIG. 3. Decay activities of isotopes within the 
228 Th decay chain. The initial build-up of daughter 
isotopes is visible in the top plot, while the bottom 
plot shows the long-term behavior of the decays. 
The vertical dashed line marks 10 times the half 
life of 224 Ra. For reference, the half life of 224 Ra is 
3.6 days. 

B. Using a 228 Th source 

Figure [3] demonstrates the validity of our approach when the top-most isotope's half life 
is short relative to the time window. The daughter isotopes are clearly populated from the 
starting parent population. The exponential trend in the overall decay rate, dictated by the 
1.9-year half life of 228 Th, is also clear. For these calculations, a starting decay activity of 
0.1 Bq was used, and 50 million decays were recorded. 
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i C. Benchmarking 



2 We benchmarked the performance of our decay generator. Table [I] shows the time required 

3 to perform the requested calculations. The calculations were run on a desktop with a 2.26 

4 GHz Intel E5520 Xeon chip and 8 GB of RAM, running MacOSX 10.5.8. The vast majority 

5 of the memory usage was taken up by the binary search tree data structure, where each node 
e of which required 30 bytes of RAM. 20 levels of pre-seeded nodes required approximately 1 
7 million entries, and the rest of the entries were specified at run time. Thus a data record 
s of 50 million events required approximately 1.5 GB of memory. This is a large amount, 

9 certainly, but still within the capabilities of any modern computer. 

10 Accurate extrapolations of the computing time required to create the event record are 

11 not straightforward, as they involve a non-trivial interplay between the pre-seeded binary 

12 search tree, the rate of primary decays, and the required number of traversals of a full decay 

13 chain. Broadly speaking, however, the computing time scales with the number of required 
u events, and the longest computing time occurs when the source age is 0. 

is IV. CONCLUSIONS 

16 The approach described in this paper for decay chain generation can be used to calculate 

17 a full, time-ordered history of decays from a radioactive source with an extended decay 
is chain, while preserving both time and position correlations between decays within a single 

19 traversal of the chain. Once these calculations are performed, the data record can be passed 

20 to a Monte Carlo program to produce the decay products themselves at the recorded times 

21 and positions. Such an event generator is useful for both basic science simulations as well 

22 as detector component evaluation, and does not make any assumptions regarding detector 

23 response time or data acquisition performance. As such, it is appropriate for a wide variety 

24 of applications. 

25 Using a purely stochastic or a purely analytic approach to decay chain generation both 

26 contain strong disadvantages. The two approaches can be combined, however, to great 

27 effect, with the analytics used to generate starting populations among the decay nuclei, and 

28 stochastics used to propagate temporal and spatial correlations. Using the combination, 

29 a large number of time-ordered decays can be generated in just a few minutes, with the 
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TABLE I. Benchmarking the decay chain generator. The second row in the headings contains the 
number of decays in the final data record. The third row in the headings is the source age, where 
the value of Ty 2 is 5.75 years, and 3.6 days for Ty 2 . All entries below the headings are in seconds, 
except for the column on the left. The first row of entries for the 228 Th decays is absent because 
the starting activity implies a starting population too low to provide the full set of decays. Each 
entry in the table represents the average and standard deviation from 10 trials. See text for details. 
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19.9 ±0.9 
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N/A 
93.8 ± 1.1 
520 ± 15 


N/A 
91.9 ± 1.3 
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computational overhead being driven primarily by the requested number of decays in the 
record. 

The data structure of the record must be chosen to minimize both search and inser- 
tion/deletion operation overhead, and thus the authors recommend a binary search tree. 
To minimize degeneracy, the search tree can be pre-seeded with empty nodes to apply a 
structure within which the decays can be recorded. 

The authors would like to thank Jason Detwiler for reviewing the manuscript, and Aric 
Stewart for recommending binary search trees for this application. Prepared by LLNL under 
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